-- ===========================
-- Minimum enclosing circle algorithm
-- ---------------------------
-- Much of this code was inspired from the ruby implementation at:
-- [http://www.dseifert.net/code/mec/source.html]
-- ===========================

local debug = true

Point = {
	-- ===========================
	-- Sets the metatable of the "point" table to
	-- the Point table. Kind of like dynamic inheritance.
	-- ===========================
	New = function(self, point)
		-- Make sure we have an actual point.
		if not point or type(point) ~= "table" or point.x == nil or point.y == nil then
			print("Point:New - point is not valid.")
			return
		end
		setmetatable(point, {__index = self, 
								__add = self.Add, 
								__sub = self.Sub,
								__eq = self.Equals,
								__tostring = self.ToString})
		return point
	end,
	
	-- ===========================
	-- Overload the "+" operator
	-- ===========================
	Add = function(self, other) 
		return Point:New({x = self.x + other.x, y = self.y + other.y}) 
	end,
	
	-- ===========================
	-- Overload the "-" operator
	-- ===========================
	Sub = function(self, other)
		return Point:New({x = self.x - other.x, y = self.y - other.y})
	end,
	
	-- ===========================
	-- Check if the other point is the same as self.
	-- ===========================
	Equals = function(self, other)
		return self.x == other.x and self.y == other.y
	end,

	-- ===========================
	-- [Static] Check to see if x and y are close, as defined by eps (if it is passed)
	-- ===========================
	CloseTo = function(x, y, eps)
		if not eps then 
			eps = 1e-9 
		end
		return math.abs(x - y) <= eps
	end,
	
	-- ===========================
	-- [Static] Get the intersection point of lines L1=(base1,v1) L2=(base2,v2)
	-- ===========================
	InfLineIntersection = function(base1, v1, base2, v2)
		if Point.CloseTo(v1.x, 0.0) and Point.CloseTo(v2.y, 0.0) then
			return Point:New({x = base1.x, y = base2.y})
		end
		if Point.CloseTo(v1.y, 0.0) and Point.CloseTo(v2.x, 0.0) then
			return Point:New({x = base2.x, y = base1.y})
		end
		
		local m1 = (not Point.CloseTo(v1.x, 0.0)) and v1.y / v1.x or 0.0
		local m2 = (not Point.CloseTo(v2.x, 0.0)) and v2.y / v2.x or 0.0
		
		if Point.CloseTo(m1, m2) then
			return nil
		end
		
		local c1 = base1.y - m1 * base1.x
		local c2 = base2.y - m2 * base2.x
		
		local ix = (c2 - c1) / (m1 - m2)
		local iy = m1 * ix + c1
		if Point.CloseTo(v1.x, 0.0) then
			return Point:New({x = base1.x, y = base1.x * m2 + c2})
		end
		if Point.CloseTo(v2.x, 0.0) then
			return Point:New({x = base2.x, y = base2.x * m1 + c1})
		end
		return Point:New({x = ix, y = iy})
	end,
	
	-- ===========================
	-- Returns a positive number if self > other, 0 if self == other, and
	-- a negative number if self < other.
	-- ===========================
	CompareTo = function(self, other)
		local ret = self.x - other.x
		if ret == 0 then
			ret = self.y - other.y
		end
		return ret
	end,
	
	-- ===========================
	-- Scale the current point to get a new point.
	-- ===========================
	Scale = function(self, factor)
		return Point:New({x = self.x * factor, y = self.y * factor})
	end,

	-- ===========================
	-- Get the string representation of the point.
	-- ===========================	
	ToString = function(self)
		return "(" .. tostring(self.x) .. ", " .. tostring(self.y) .. ")"
	end,
	
	-- ===========================
	-- Computes the angle formed by p1 - self - p2
	-- ===========================
	AngleBetween = function(self, p1, p2)
		local vect_p1 = p1 - self
		local vect_p2 = p2 - self
		local theta = vect_p1:Polar() - vect_p2:Polar()
		if theta < 0.0 then
			theta = theta + 360.0
		end
		if theta > 180.0 then
			theta = 360.0 - theta
		end
		return theta
	end,
	
	-- ===========================
	-- Compute the polar angle to this point
	-- ===========================
	Polar = function(self)
		theta = 0.0
		if Point.CloseTo(self.x, 0) then
			if self.y > 0 then
				theta = 90.0
			elseif self.y < 0 then
				theta = 270.0
			end
		else
			theta = math.deg(math.atan(self.y/self.x))
			if self.x < 0.0 then
				theta = theta + 180.0
			end
			if theta < 0.0 then
				theta = theta + 360.0
			end
		end
		return theta
	end,
	
	-- ===========================
	-- Gets the left normal according to self's position.
	-- ===========================
	NormalLeft = function(self)
		return Point:New({x = self.y, y = -self.x})
	end,
	
	-- ===========================
	-- Get the center between self and the other point.
	-- ===========================
	Center = function(self, otherPoint)
		return Point:New({x = (self.x + otherPoint.x) / 2, y = (self.y + otherPoint.y) / 2})
	end,
	
	-- ===========================
	-- Get the distance between myself and the other point.
	-- ===========================
	Distance = function(self, otherPoint)
		local dx = self.x - otherPoint.x
		local dy = self.y - otherPoint.y
		return math.sqrt(dx * dx + dy * dy)
	end,
	
	-- ===========================
	-- Tests if a point is left|on|right of an infinite line
	-- return >0 for p2 left of line through p0 and p1.
	-- return =0 for p2 on the line.
	-- return <0 for p2 right of the line.
	-- ===========================
	Direction = function(p0, p1, p2)
		return (p0.x-p1.x)*(p2.y-p1.y) - (p2.x-p1.x)*(p0.y-p1.y)
	end,
}

Circle = {
	-- ===========================
	-- Creates a new circle of the specified center and radius.
	-- ===========================
	New = function(self)	
		local c = {}
		setmetatable(c, {__index = self})
		return c
	end,
	
	-- ===========================
	-- Determines if the point is contained within the circle.
	-- ===========================
	Contains = function(self, point)
		return point:Distance(self.center) < self.radius 
				or math.abs(point:Distance(self.center) - self.radius) < 1e-9
	end,
	
	-- ===========================
	-- Get the string representation of the circle.
	-- ===========================
	ToString = function(self)
		return "{center: " .. self.center:ToString() .. ", radius: " .. tostring(self.radius) .. "}"
	end,
}

-- ===========================
-- Copy the array portion of the table.
-- ===========================
function CopyArray(t)
	local ret = {}
	for i,v in ipairs(t) do
		ret[i] = v
	end
	return ret
end

-- ===========================
-- Create a new array, minus the element at the specified index
-- ===========================
function DeleteAt(t, index)
	local ret = {}
	for i,v in ipairs(t) do
		if i ~= index then
			table.insert(ret, v)
		end
	end
	return ret
end

-- ===========================
-- Joins arrays t1 and t2 and outputs only the unique elements.
-- ===========================
function JoinUniqueArrayElements(t1, t2)
	local ret = {}
	local unique = {}
	-- Place the elements of t1 and t2 into the unique dictionary
	for i,p in ipairs(t1) do
		unique["x" .. p.x .. "y" .. p.y] = p 
	end
	for i,p in ipairs(t2) do
		unique["x" .. p.x .. "y" .. p.y] = p
	end
	-- Insert each element of unique into the return array.
	for k,p in pairs(unique) do
		table.insert(ret, p)
	end
	return ret
end

-- ===========================
-- Minimum Enclosing Circle object and algorithm
-- ===========================
MEC = {
	-- ===========================
	-- Create a new MEC table
	-- ===========================
	New = function(self, points)
		local mec = {}
		mec = setmetatable({}, {__index = self})
		mec.circle = nil
		mec.points = {} -- a table of x,y coordinates
		
		if points then
			mec:SetPoints(points)
		end
		
		return mec
	end,
	
	-- ===========================
	-- Set the points used to compute the MEC.
	-- points is an array, starting at 1.
	-- ===========================
	SetPoints = function(self, points)
		-- Set the points
		self.points = points
		for i, p in ipairs(self.points) do
			p = Point:New(p)
		end
	end,
	
	-- ===========================
	-- Computes the half hull of a set of points
	-- ===========================
	HalfHull = function(left, right, pointTable, factor)
		local input = CopyArray(pointTable)
		table.insert(input, right)
		local half = {}
		table.insert(half, left)
		for i,p in ipairs(input) do
			table.insert(half, p)
			while #half >= 3 do
				local dir = factor * Point.Direction(half[(#half+1)-3], half[(#half+1)-1], half[(#half+1)-2])
				if dir <= 0 then
					half = DeleteAt(half, #half-1)
				else
					break
				end
			end
		end
		return half
	end,
	
	-- ===========================
	-- Computes the set of points that represent the
	-- convex hull of the set of points
	-- ===========================
	ConvexHull = function(self)
		local a = self.points
		local left = a[1]
		local right = a[#a]
		local upper = {}
		local lower = {}
		
		-- Partition remaining points into upper and lower buckets.
		for i = 2, #a-1 do
			local dir = Point.Direction(left, right, a[i])
			if dir < 0 then
				table.insert(upper, a[i])
			else
				table.insert(lower, a[i])
			end
		end
		
		local upperHull = self.HalfHull(left, right, upper, -1)
		local lowerHull = self.HalfHull(left, right, lower, 1)
		
		return JoinUniqueArrayElements(upperHull, lowerHull)
	end,
	
	-- ===========================
	-- Compute the MEC.
	-- ===========================
	Compute = function(self)
		self.circle = self.circle or Circle:New()
		
		-- Make sure there are some points.
		if #self.points == 0 then return self.circle end
		
		-- Handle degenerate cases first
		if #self.points == 1 then
			self.circle.center = self.points[1]
			self.circle.radius = 0
			self.circle.radiusPoint = self.points[1]
		elseif #self.points == 2 then
			local a = self.points
			self.circle.center = a[1]:Center(a[2])
			self.circle.radius = a[1]:Distance(self.circle.center)
			self.circle.radiusPoint = a[1]
		else
			local a = self:ConvexHull()
			local point_a = a[1]
			local point_b = nil
			local point_c = a[2]
			
			if not point_c then
				self.circle.center = point_a
				self.circle.radius = 0
				self.circle.radiusPoint = point_a
				return self.circle
			end
			
			-- Loop until we get appropriate values for point_a and point_c
			while true do
				point_b = nil
				local best_theta = 180.0
				-- Search for the point "b" which subtends the smallest angle a-b-c. 
				for i,point in ipairs(self.points) do
					if (not point:Equals(point_a)) and (not point:Equals(point_c)) then
						local theta_abc = point:AngleBetween(point_a, point_c)
						if theta_abc < best_theta then
							point_b = point
							best_theta = theta_abc
						end				
					end
				end
				-- If the angle is obtuse, then line a-c is the diameter of the circle,
				-- so we can return.
				if best_theta >= 90.0 or (not point_b) then
					self.circle.center = point_a:Center(point_c)
					self.circle.radius = point_a:Distance(self.circle.center)
					self.circle.radiusPoint = point_a
					return self.circle
				end
				local ang_bca = point_c:AngleBetween(point_b, point_a)
				local ang_cab = point_a:AngleBetween(point_c, point_b)
				if ang_bca > 90.0 then
					point_c = point_b
				elseif ang_cab <= 90.0 then
					break
				else
					point_a = point_b
				end
			end
			local ch1 = (point_b - point_a):Scale(0.5)
			local ch2 = (point_c - point_a):Scale(0.5)
			local n1 = ch1:NormalLeft()
			local n2 = ch2:NormalLeft()
			ch1 = point_a + ch1
			ch2 = point_a + ch2
			self.circle.center = Point.InfLineIntersection (ch1, n1, ch2, n2)
			self.circle.radius = self.circle.center:Distance(point_a)
			self.circle.radiusPoint = point_a
		end
		return self.circle
	end,
}


